home *** CD-ROM | disk | FTP | other *** search
/ NetNews Offline 2 / NetNews Offline Volume 2.iso / news / comp / lang / c-part1 / 3327 < prev    next >
Encoding:
Text File  |  1996-08-05  |  2.3 KB  |  92 lines
  1. Path: ocbbs.gen.nz!not-for-mail
  2. From: steve@hn.ocbbs.gen.nz (Steve Detoni)
  3. Newsgroups: comp.lang.c
  4. Subject: Re: Fast Integer Square Root Algorithm
  5. Date: 28 Jan 1996 11:11:00 +1300
  6. Message-ID: <4ee7tk$7lk@hn.ocbbs.gen.nz>
  7. References: <cpruett-2001962141520001@osip74.ionet.net>
  8. NNTP-Posting-Host: hn.hn.planet.gen.nz
  9. X-Newsreader: TIN [version 1.2 PL2]
  10.  
  11. Chris Pruett (cpruett@ionet.net) wrote:
  12. : Can anyone lend suggestions or give me a reference to a very
  13. : fast integer square root algorithm?
  14.  
  15. : 1) Must take square root of unsigned, long integer (32-bit).
  16. : 2) Accuracy is not terribly important.  +/- 100 is acceptable.
  17.  
  18. : The best I've been able to come up with is Newton's method.  
  19. : Is there something better?  
  20. I asked this question myself, and got a reply, the following code was 
  21. posted me and i'm sure its the fastest method. Doesn't use multiplcation, 
  22. but successive approimatation using bit-shifting. Tis very fast.
  23.  
  24. #include <stdio.h>
  25. #include <stdlib.h>
  26.  
  27. #define BITS_PER_BYTE 8      // Just in case, but it better be even
  28. #define MYTYPE unsigned long // For generality?
  29.  
  30.  
  31. void main( int argc, char *argv[])
  32. {
  33.   MYTYPE question;
  34.   MYTYPE mask;
  35.   MYTYPE answer;
  36.   MYTYPE temp;
  37.   MYTYPE divisor;
  38.   int rotval;
  39.  
  40.   // Just for the purpose of demo
  41.   if (argc == 2)
  42.   {
  43.         question = (MYTYPE)atol( argv[ 1]); // Fudge here for demo
  44.   }
  45.   else
  46.   {
  47.         printf(" Usage %s <number>\n", argv[ 0]);
  48.         exit( 1);
  49.   }
  50.  
  51.   answer  = (MYTYPE)0;
  52.   divisor = (MYTYPE)0;
  53.   rotval  = ((sizeof(question) * BITS_PER_BYTE) - 2);
  54.   mask    = (MYTYPE)3 << rotval;
  55.   temp    = (MYTYPE)0;
  56.  
  57.   temp    = (temp << 2) + ((question & mask) >> rotval);
  58.   answer  = (answer << 1) + (divisor & 1);
  59.   divisor = ((divisor + (divisor & 1))<<1) + 1;
  60.  
  61.   do
  62.   {
  63.       if (divisor > temp)
  64.       {
  65.           divisor--;
  66.       }
  67.       else
  68.       {
  69.           temp -= divisor;
  70.       }
  71.  
  72.       rotval -= 2;
  73.       mask = mask >> 2;
  74.  
  75.       temp    = (temp << 2)   + ((question & mask) >> rotval);
  76.       answer  = (answer << 1) + (divisor & 1);
  77.       divisor = ((divisor     + (divisor & 1))<<1) + 1;
  78.  
  79.   } while (mask);
  80.  
  81.   // Round up if required (next bit would be .1 ie >= .5 in decimal
  82.   if (temp >= divisor)
  83.   {
  84.       answer++;
  85.   }
  86.  
  87.   printf("Answer = %ld\n", answer);     // Not very general :)
  88. } 
  89.  
  90. Hope this helps.
  91. Steve.
  92.